A Measure of Nonlinear Correlation with an Estimation Algorithm
従属性の指標が満たすべき性質はTasena & Dhompongsa (2016)などで
議論されているが、当面は素朴な定義を用いる
連続変数の値域を分割し、離散変数として評価する
\[ f(x,y) \]
\[ p(i,j) = P(X \in \mathcal{X}_i, Y \in \mathcal{Y}_j) \]
\[ X \in \mathcal{X}, \; Y \in \mathcal{Y} \]
\[ \mathcal{X} = \sqcup_{i=1}^k \mathcal{X}_i, \; \mathcal{Y} = \sqcup_{j=1}^k \mathcal{Y}_j \]
\[ \begin{align*} I(X;Y) &= \mathbb{E} \Bigg[ \log \frac{p(x,y)}{p_X(x)p_Y(y)} \Bigg] \\ &= \sum_x \sum_y p(x,y) \log \frac{p(x,y)}{p_X(x)p_Y(y)} \end{align*} \]
\[ 0 \le I(X;Y) \le H(Y) \]
\[ H(Y) = - \sum_y p_Y(y) \log p_Y(y) \]
相互情報量を目的変数の情報量で基準化
\[ U(Y|X) = \frac{I(X;Y)}{H(Y)} \]
\[ 0 \le U(Y|X) \le 1 \]
\(X\)によって説明できる、\(Y\)の不確実性の割合
探索的データ解析ソフトウェアNattoにおいて、変数間の関係を評価する手法として実装 (Suzuki et al., 2006) https://ef-prime.com/ja/product/
上記の例は相互情報量や不確実性係数の意味で等価
\[ \begin{align*} I(X;Y) &= \mathbb{E} \Bigg[ \log \frac{f(x,y)}{f_X(x)f_Y(y)} \Bigg] \\ &= \iint f(x,y) \log \frac{f(x,y)}{f_X(x)f_Y(y)} dx dy \end{align*} \]
以下、\(I\)と略記する
\(X, Y\)の値域をそれぞれ\(k\)分割した確率変数\(X_k, Y_k\)の相互情報量
\[ \begin{align*} I(X_k;Y_k) &= \mathbb{E} \Bigg[ \log \frac{p(i,j)}{p_X(i)p_Y(j)} \Bigg] \\ &= \sum_i^k \sum_j^k p(i,j) \log \frac{p(i,j)}{p_X(i)p_Y(j)} \end{align*} \]
これを\(I_k\)と書くと\(I_k \le I\)
\(I\)の被積分関数がリーマン可積分なら、分割を細かくすれば\(I_k \to I \; (k \to \infty)\)
不確実係数\(U(Y_k|X_k)\)について、\(k \to \infty\)の挙動を考える:
\[ U(Y_k|X_k) = \frac{I(X_k;Y_k)}{H(Y_k)} \]
\[ 0 \le U(Y_k|X_k) \le 1 \]
MIC (Reshef et al., 2011) も変数を離散化し、相互情報量を基準化する。基準化した値の最大値をとることで、分母の発散に対応していると解釈できる \[ \max_{s, t} \frac{\hat{I}(X_s;Y_t)}{\log (\min\{s, t\})} \]
理論値としての評価指標を定め、その推定手法を考える
\[ I(X;Y) = \begin{cases} \displaystyle \mathbb{E}_{P_{X,Y}} \Bigg[ \log \frac{dP_{X,Y}}{dP_X \otimes dP_Y} \Bigg] & (P_{X,Y} \ll P_X \otimes P_Y) \\ \infty & (\text{otherwise}) \end{cases} \]
\(dP/dQ\)はラドン=ニコディム微分
対数変換の扱いに工夫を要する場合がある。例として推定量の分布を
評価する際に多項式近似が用いられる (Hamdan & Tsokos, 1971)
\[ \psi(X;Y) := \begin{cases} \displaystyle \mathbb{E}_{P_{X,Y}} \Bigg[ \frac{dP_{X,Y}}{dP_X \otimes dP_Y} \Bigg] & (P_{X,Y} \ll P_X \otimes P_Y) \\ \infty & (\text{otherwise}) \end{cases} \]
相互情報量から対数変換を除いたものとして定義
逆数は\(\psi^{-1} \in [0,1]\)で、独立性の指標となり得る
対数変換を含まないことで、式展開や推定量の評価が容易になる
\[ \psi(X;Y) = \mathbb{E}_{P_{X,Y}} \Bigg[ \frac{f(x,y)}{f_X(x)f_Y(y)} \Bigg] = \mathbb{E}_{P_{X,Y}} \Bigg[ \frac{f_{Y|X}(y|x)}{f_Y(y)} \Bigg] \]
Jensenの不等式による相互情報量の上界を構成する:
\[ \begin{align*} I(X;Y) &= \mathbb{E}_{P_{X,Y}} \Bigg[ \log \frac{dP_{X,Y}}{dP_X \otimes dP_Y} \Bigg]\\ &\le \log \Bigg( \mathbb{E}_{P_{X,Y}} \Bigg[ \frac{dP_{X,Y}}{dP_X \otimes dP_Y} \Bigg] \Bigg) \\ &= \log \psi(X;Y) \end{align*} \]
\[ \begin{alignat*}{1} I(X;Y) &= D_{\rm KL} & (P_{X,Y} \Vert P_X \otimes P_Y) \\ \psi(X;Y) &= D_{\chi^2} & (P_{X,Y} \Vert P_X \otimes P_Y) + 1 \end{alignat*} \]
\(\Rightarrow\) \(f\)-ダイバージェンスや\(\chi^2\)統計量の性質が利用できる
相互情報量と同様の望ましい性質をもつ:
相互情報量を用いても類似の議論が可能だが、
対数変換を含まないことで各種の議論がシンプルになる
\((X,Y)\)が同時密度\(f(x,y)\)、周辺密度\(f_X(x), f_Y(y)\)をもつとき
\[ \begin{align*} \psi(X;Y) &= \mathbb{E}_{P_{X,Y}} \Bigg[ \frac{f(x,y)}{f_X(x)f_Y(y)} \Bigg] \\ &= \iint \frac{f(x,y)^2}{f_X(x)f_Y(y)} dx dy \end{align*} \]
\((X,Y)\)が相関係数\(\rho\)をもつ2変量正規分布に従うとき
\[ \psi = \frac{1}{1-\rho^2} \]
が成り立ち、\(\psi\)と\(\rho^2\)は一対一に対応する。すなわち
\[ \rho^2 = 1 - \psi^{-1} \]
正規分布における議論に基づいて以下を定義(*):
\[ r_\psi = \sqrt{1 - \psi^{-1}} \]
* 相互情報量に基づく Linfoot (1957) の情報相関係数
(informational coefficient of correlation)と同様の発想\[ r_I = \sqrt{1 - e^{-2 I(X;Y)}} \]
\((X,Y)\)が同時確率\(p(x,y)\)、周辺確率\(p_X(x), p_Y(y)\)をもつとき
\[ \begin{align*} \psi(X;Y) &= \mathbb{E}_{P_{X,Y}} \Bigg[ \frac{p(x,y)}{p_X(x)p_Y(y)} \Bigg] \\ &= \sum_x \sum_y \frac{p(x,y)^2}{p_X(x)p_Y(y)} \end{align*} \]
\(X,Y\)の支持集合の基数をそれぞれ\(k_X,k_Y\)とするとき
\[ 1 \le \psi(X;Y) \le \min \{k_X, k_Y\} \]
\(Y\)が\(X\)に完全従属する、すなわち\(p_{Y|X}\)が一点分布となるとき
\[ \psi(X;Y) = k_Y \le k_X \]
一方が連続または無限離散などの場合も同様の議論が可能
相互依存度に基づく予測スコア
(Predictability Score)を以下で定義する:
\[ \mathrm{Pred}_\psi(X \to Y) := \begin{cases} \displaystyle \sqrt{\frac{1-\psi^{-1}}{1-k_Y^{-1}}} & (k_Y > 1) \\[1em] 1 & (k_Y = 1) \end{cases} \]
\(k_Y\)は\(Y\)の支持集合の基数(無限集合の場合は\(\infty\))
予測スコアは予測可能性の非対称性が反映されているが、
相関係数のように対称性をもつ指標があると便利
以下を満たす指標を一般化相関尺度(Generalized Correlation Measures)
と呼ぶことにする:
相互依存度に基づく一般化相関尺度を以下で定義する:
\[ \mathrm{gCor}_\psi(X,Y) := \begin{cases} \displaystyle \sqrt{\frac{1-\psi^{-1}}{1- \big(\sqrt{k_Xk_Y}\big)^{-1}}} &(k_X k_Y > 1) \\[1em] 1 & (k_X k_Y = 1) \end{cases} \]
\(k_X, k_Y\)は\(X, Y\)の支持集合の基数(無限集合の場合は\(\infty\))
以降\(\mathrm{gCor}(X,Y)\)または単に\(\mathrm{gCor}\)と略記
連続確率変数\(X, Y\)の値域をそれぞれ\(k\)分割した\(X_k, Y_k\)の相互依存度
\[ \begin{align*} \psi(X_k;Y_k) &= \mathbb{E} \Bigg[ \frac{p(i,j)}{p_X(i)p_Y(j)} \Bigg] \\ &= \sum_i^k \sum_j^k \frac{p(i,j)^2}{p_X(i)p_Y(j)} \end{align*} \]
これを\(\psi_k\)と書くと\(\psi_k \le \psi\)
\(\psi\)の被積分関数がリーマン可積分なら\(\psi_k \to \psi \; (k \to \infty)\)
一般化相関尺度(および予測スコア)についても同様:
\[ \mathrm{gCor}(X_k,Y_k) \to \mathrm{gCor}(X,Y) \;\; (k \to \infty) \]
相関係数\(\rho = 0.5\)をもつ2変量正規分布を\(k\)分位点で分割
○が離散化による近似値\(\psi_k\)、赤線が真値\(\psi = 3/4\)
○は離散化による近似値\(\psi_k\)を式\(\sqrt{1 - \psi^{-1}}\)に代入した値
×が\(\sqrt{1 - k^{-1}}\)で基準化した\(\mathrm{gCor}\)の値、赤線が真値\(|\rho| = 0.5\)
離散確率変数\(V,W\)が(観測されない)連続確率変数\(X,Y\)の分割
\[ V = v(X), \ W = w(Y) \]
と仮定すれば、\(\mathrm{gCor}(V,W)\)を\(\mathrm{gCor}(X,Y)\)の近似として解釈できる
特に\((X,Y)\)が2変量正規分布に従うとき、相関係数の絶対値\(|\rho|\)に対する近似となる
\(\Rightarrow\)連続、離散および混合型に対する統一的な評価に意味を与える
例:\(\mathrm{gCor}(V,W) = 0.7\)のとき、背後に\(\mathrm{Cor}(X,Y) = 0.7\)となる正規確率変数\(X,Y\)がある状況に近いと考える
数値データ(実数または整数)を分割して離散化
離散化データに対して\(\mathrm{gCor}\)などの指標を推定
プラグイン推定量\(\hat\psi\)は\(\chi^2\)統計量で表すことができる:
\[ \hat{\psi} = \sum_i \sum_j \frac{n^2_{ij}}{n_{i.}n_{.j}} = \frac{\chi^2}{n} + 1 \]
\(\Rightarrow\) 近似的に\(\chi^2\)分布を用いた統計的推測が可能
\[ n_{ij} = \#(X = x_i, Y = y_j), \ \ n_{i.} = \#(X = x_i), \ \ n_{.j} = \#(Y = y_j)\\ n = \sum_i \sum_j n_{ij}, \ \ \chi^2 = \sum_i \sum_j \frac{(n_{ij} - n_{i.}n_{.j}/n)^2}{n_{i.}n_{.j}/n} \]
\(X,Y\)が離散のとき、クラメールの\(V\)は\(\chi^2\)統計量で表される:
\[ V = \sqrt{\frac{\chi^2/n}{\min\{k_X, k_Y\} - 1}} \]
\(k_X = k_Y = 2\)のときピアソンの\(\phi\)が定義され、\(V = |\phi|\)となる
相互依存度のプラグイン推定量は\(V\)の分子に1を足した値となり、
逆数をとって変換することで一般化相関尺度の推定値が得られる:
\[ \hat{\psi} = \frac{\chi^2}{n} + 1 \]
\[ \widehat{\mathrm{gCor}}(X,Y) = \sqrt{\frac{1-\hat{\psi}^{-1}}{1-\big(\sqrt{k_Xk_Y}\big)^{-1}}} \]
\((X,Y)\)が相関係数\(\rho\)をもつ2変量正規分布に従うとき、\(k\)分割による離散化を行えば、適当な条件のもとで以下が得られる:
\[ \begin{align*} \sqrt{1 - \hat\psi_k^{-1}} &= \sqrt{1 - \bigg(\frac{\chi^2}{n} + 1\bigg)^{-1}} \\[1em] &= \sqrt{\frac{\chi^2/n}{1 + \chi^2/n}} \\[1em] &\to |\rho| \ \ \ (k,n \to \infty) \end{align*} \]
K. Pearson (1904) において\(\phi^2 = \chi^2/n\)をmean square contingencyと定義し、類似の議論を展開している。上記に相当する値をfirst coefficient of contingencyとして提案しており、\(\widehat{\mathrm{gCor}}\)はこの値を分割数で
基準化したものにあたる
\[ \max \bigg\{ 2, \big\lfloor \sqrt{n / 50} \big\rfloor \bigg\} \]
以下を1,000回繰り返す:
\(\Rightarrow\) 推定値の平均値、および90%区間をプロットして確認
以降、セルあたり平均サンプルサイズを平均サイズと記載
\(\sqrt{1000/50} \approx 4.47\)を参考に\(k\)を設定、平均サイズ62.5
赤の点線は真値。バイアスはあるものの真値に近く、バラツキも小さい
\(\sqrt{100/50} \approx 1.41\)を参考に\(k\)を設定、平均サイズ25
バラツキが大きくなるものの、バイアスはそこまで増加しない
\(\sqrt{10000/50} \approx 14.14\)を参考に\(k\)を設定、平均サイズ約51
バイアスは残るが、バラツキはほとんどなくなる
\(\sqrt{100/50} \approx 1.41\)を大幅に上回り、平均サイズは6.25
バイアスは大きくなるが、平均的な順序関係は保たれている
平均サイズが1となる極端な設定
バイアスが非常に大きくなり、バラツキはかえって小さくなる
\(\sqrt{1000/50} \approx 4.47\)に対して約2.2倍に\(k\)を設定、平均サイズ10
バイアスの傾向はn = 100, k = 4の例(平均サイズ6.25)と類似
\(n\)を10倍に増やし、\(\sqrt{10000/50} \approx 14.14\)に対して約0.7倍に\(k\)を設定
平均サイズは100に増加し、バイアスが大きく改善
\(k\)を2倍に増加し、平均サイズは25に減少
バイアスが増加するが、非線形な関係を検知しやすくなると期待
非線形な関係を検出するためには分割を増やす必要がある
バイアスを許容して分割を増やし、関係の検知性能を重視する方針も
ランダム生成した2次元の数値データを評価
サンプルサイズ\(n\)、分割数\(k\)を動かして挙動を確認
標本相関係数と比較
Wikipedia「ピアソンの積率相関係数」掲載の例をもとに作成
上段は相関係数\(\rho\)の絶対値を1, 0.8, 0.4, 0に設定した2変量正規分布
小数点以下は1桁に丸めて表示(以降も同様)
平均サイズ50を基準として\(k\)を設定(\(n/k^2 = 62.5\))
非線形なパターンを検出しつつ、正規分布に対しては\(|\rho|\)に近い値
サンプルサイズが小さい場合
正規乱数に対しても多少の推定誤差が生じる
平均サイズ50を基準として\(k\)を設定(\(n/k^2 = 25\))
正規分布に対しては\(|\rho|\)に近いが、非線形パターンを検出できていない
平均サイズ50を大幅に下回る(\(n/k^2 \approx 11\))
近似バイアスが増加するものの、非線形パターンの検出には成功
サンプルサイズが大きい場合について、小数点以下第2位まで表示
サンプルサイズを増やしても、非線形パターンについては0に近い
平均サイズ50を基準として\(k\)を設定(\(n/k^2 \approx 51\))
近似バイアスが確認できるものの、非線形パターンの検出は良好
平均サイズ100を基準として\(k\)を設定(\(n/k^2 = 100\))
独立に近い例で近似バイアスが低下している
一方で以下の点は提案手法のメリットといえる:
離散化された変数を\(X_s, Y_t\)とおき、支持集合の基数を\(s,t\)とする
\(\chi^2\)統計量を\(\hat{\chi}^2\)と書くと、独立性の帰無仮説のもとで
\[ \hat{\chi}^2 \sim \chi^2_{(s - 1)(t - 1)} \; \; (n \to \infty) \]
このとき期待値は\(E\big[\hat{\chi}^2\big] \approx (s - 1)(t - 1)\)
プラグイン推定量\(\hat{\psi} = \hat{\chi}^2 / n + 1\)について
\[ E \big[ \hat{\psi} \big] \approx 1 + \frac{(s - 1)(t - 1)}{n} \]
帰無仮説のもとで\(\psi = 1\)
\(\Rightarrow\) 独立であっても分割数が増えると上方バイアスが生じる
帰無仮説のもとで(漸近的に)不偏な推定量を定義:
\[ \tilde{\psi} \; := \; \hat{\psi} - \frac{(s - 1)(t - 1)}{n} \]
\(\psi \ge 1\)より、推定値としては\(\max\{1, \tilde{\psi}\}\)を用いる
- クラメールのVに関しても類似のバイアス補正が提案されている (Bergsma, 2013)
- 正確確率検定(周辺度数を固定)のもとで分母は\(n-1\)
gcor開発版をGitHubにて公開中(オープンソース)
https://github.com/r-suzuki/gcor
irisデータ
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1 5.1 3.5 1.4 0.2 setosa
2 4.9 3.0 1.4 0.2 setosa
3 4.7 3.2 1.3 0.2 setosa
4 4.6 3.1 1.5 0.2 setosa
5 5.0 3.6 1.4 0.2 setosa
6 5.4 3.9 1.7 0.4 setosa
情報をもつ欠損(informative missingness)
Species == "setosa"のとき50%の確率でSepal.Widthが欠損
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1 5.1 NA 1.4 0.2 setosa
2 4.9 3.0 1.4 0.2 setosa
3 4.7 NA 1.3 0.2 setosa
4 4.6 3.1 1.5 0.2 setosa
5 5.0 3.6 1.4 0.2 setosa
6 5.4 NA 1.7 0.4 setosa
カテゴリ値を含むためエラー
数値変数のみに限定して実行
数値とカテゴリを統一的に評価
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
Sepal.Length 1.0000000 0.2349075 0.8846517 0.8741873 0.7623968
Sepal.Width 0.2349075 1.0000000 0.3143301 0.2669031 0.6510740
Petal.Length 0.8846517 0.3143301 1.0000000 0.9503289 0.8221674
Petal.Width 0.8741873 0.2669031 0.9503289 1.0000000 0.8237429
Species 0.7623968 0.6510740 0.8221674 0.8237429 1.0000000
分割数\(k\)を指定(既定ではサンプルサイズから自動設定)
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
Sepal.Length 1.0000000 0.6180293 0.8449216 0.8300394 0.8226041
Sepal.Width 0.6180293 1.0000000 0.6862189 0.6963419 0.6794982
Petal.Length 0.8449216 0.6862189 1.0000000 0.9581278 0.9728947
Petal.Width 0.8300394 0.6963419 0.9581278 1.0000000 0.9795837
Species 0.8226041 0.6794982 0.9728947 0.9795837 1.0000000
欠損を含む列は算出されない
列の組ごとに、欠損を含むケースを無視する
Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length 1.0000000000 0.0002584192 0.8717538 0.8179411
Sepal.Width 0.0002584192 1.0000000000 -0.3163403 -0.2364441
Petal.Length 0.8717537759 -0.3163403009 1.0000000 0.9628654
Petal.Width 0.8179411263 -0.2364441194 0.9628654 1.0000000
「欠損という事象を観測した」として評価
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
Sepal.Length 1.0000000 0.5017295 0.8846517 0.8741873 0.7623968
Sepal.Width 0.5017295 1.0000000 0.5347735 0.5141064 0.7313361
Petal.Length 0.8846517 0.5347735 1.0000000 0.9503289 0.8221674
Petal.Width 0.8741873 0.5141064 0.9503289 1.0000000 0.8237429
Species 0.7623968 0.7313361 0.8221674 0.8237429 1.0000000
欠損を無視することも可能(この例では評価が変わる)
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
Sepal.Length 1.000000000 0.009496501 0.88465174 0.87418733 0.7623968
Sepal.Width 0.009496501 1.000000000 0.09600307 0.04385796 0.5701664
Petal.Length 0.884651737 0.096003072 1.00000000 0.95032889 0.8221674
Petal.Width 0.874187331 0.043857963 0.95032889 1.00000000 0.8237429
Species 0.762396795 0.570166409 0.82216737 0.82374290 1.0000000
\(\Rightarrow\) Rによるソフトウェア実装を公開し、容易に実行可能
\(\Rightarrow\) いずれの観点からも近似バイアスが課題